Computing topological invariants without inversion symmetry 
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We consider the problem of calculating the weak and strong topological indices in noncentrosym- 
metric time- reversal (T) invariant insulators. In 2D we use a gauge corresponding to hybrid Wannier 
functions that are maximally localized in one dimension. Although this gauge is not smoothly de- 
fined on the two-torus, it respects the T symmetry of the system and allows for a definition of the 
Z2 invariant in terms of time-reversal polarization. In 3D we apply the 2D approach to T-invariant 
planes. We illustrate the method with first-principles calculations on GeTe and on HgTe under [001] 
and [111] strain. Our approach differs from ones used previously for noncentrosymmetric materials 
and should be easier to implement in ab initio code packages. 
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I. INTRODUCTION 

A series of theoretical developments starting in 2005, 
showing that non-magnetic insulators admit a topo- 
logical Z2 classification in two dimensions (2D)^'^ and 
then in three dimensions (SD),'^'''^ has sparked enor- 
mous interest, especially after numerous realizations of 
such systems were confirmed both thcoretically^^^ and 
experimentally. ^'^^^^ These developments, nicely summa- 
rized in some recent reviews, have essentially given 
rise to a new subfield of condensed-matter physics, with 
the topology of the band structure now regarded as a 
fundamental characteristic of the electronic ground state 
for semiconductors and insulators. 

The Z2 classification divides time-reversal (7~) invari- 
ant band insulators into two classes: ordinary (Z2-even) 
insulators that can be adiabatically converted to the vac- 
uum (or to each other) without a bulk gap closure, and 
"topological" (Z2-odd) ones that cannot be so connected 
(although they can be adiabatically connected to each 
other). Even and odd phases are separated by a topo- 
logical phase transition, and the bulk gap has to van- 
ish at the transition point, at least in a non-interacting 
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The Z2-odd states are characterized by the 



presence of an odd number of Kramers pairs of counter- 
propagating edge states in 2D, or by an odd number of 
Fermi loops enclosing certain high-symmetry points of 
the surface band structure in 3D. 

In view of all this, there is an obvious motivation to 
develop simple yet effective methods for computing the 
topological indices of a given material. For centrosym- 
mctric crystals, a convenient method was introduced in 
Rcf. 6, where it was shown that the knowledge of the 
parity eigenvalues of the electronic states at only four T- 
invariant momenta in 2D (or eight of them in 3D) is suffi- 
cient to compute the topological characteristics of a given 
material. This approach is limited to centrosymmetric 
systems, however, and the calculation of the Z2 invariant 
for noncentrosymmetric insulators is not so trivial. 

One possible approach, suggested in Ref. 20, is based 
on the existence of a topological obstruction to choos- 
ing a smooth gauge that respects the T symmetry in the 



Z2-odd case. For the implementation of this method, a 
gauge must be chosen on the boundary of half of the Bril- 
louin zone (BZ) in such a way as to respect T symmetry, 
which involves acting with the time-reversal operator on 
one of the states from each Kramers pair to construct the 
other. Although this method has been implemented in 
the ab initio framework^^"^'^, its implementation is basis- 
set dependent and involves the application of a unitary 
rotation to the computed eigenvectors when fixing the 
gauge, which may be tedious when there are many occu- 
pied bands and basis states. 

Another existing method^ relies on the fact that the 
system will necessarily be in the Z2-even (normal) state 
in the absence of spin-orbit (SO) coupling. In this 
method, the strength of the SO coupling is artificially 
tuned from Xso = (no SO coupling) to Xso = 1 (full 
SO coupling), and a closure of the band gap at some in- 
termediate coupling strength is taken as evidence of an 
inverted band structure. However, a closure of the band 
gap in the course of tuning Xso to full strength is a neces- 
sary, but not a sufficient, condition for a topological phase 
transition. Therefore, in order to determine whether the 
system is really in the topologically nontrivial phase, a 
first-principles calculation of the surface states is carried 
out in order to count the number of Dirac cones at the 
surface of the candidate material. Such a calculation, al- 
though illustrative, is quite demanding in terms of com- 
putational resources. 

In summary, existing methods have some shortcom- 
ings, and it would be very useful to develop a simple 
and effective method that would use the electronic wave- 
functions, as obtained directly from the diagonalization 
procedure, to determine the desired topological indices. 

In this paper we develop a method for computing Z2 
invariants that meets these criteria, and which is easy 
to implement in the context of ab initio code packages. 
The method is based on the concept of time-reversal 
polarization^^ (TRP), but implemented in such a way 
that a visual inspection of plotted curves is not required 
in order to obtain the topological indices. Instead, all 
the indices can be obtained directly as a result of an au- 
tomated calculation. We describe the method, and then 
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verify it using centrosymmetric Bi and Bi2Se3 as illustra- 
tive test examples before applying it to the more difficult 
cases of noncentrosymmetric GeTe and strained HgTc. 

The paper is organized as follows. In Sec. II we start 
by reviewing the formalism of TRP in the context of the 
Z2 spin pump in one dimension (ID), emphasizing its re- 
lation to the charge centers of Wannicr functions. Wc 
then discuss the numerical implementation of these ideas 
to 2D and 3D cases, and suggest a simple numerical pro- 
cedure for calculating the Z2 invariant in noncentrosym- 
metric T-invariant systems in Sec. III. Wc further illus- 
trate this method with ab initio calculations in Sec. IV, 
and present some concluding remarks in Sec. V. 



II. Z2 INVARIANT VIA WANNIER CHARGE 
CENTERS 

In this section we review the notion of TRP and the 
definition of the Z2 invariant in terms of TRP derived 
in Ref. 24. The definition arises by virtue of an analogy 
between a 2D T-invariant insulator and a 7~-symmctric 
pumping process in a ID insulator. Wc further reformu- 
late this definition in terms of Wannier charge centers, 
setting the stage for the numerical method discussed in 
the next section. 



A. Review of time reversal polarization 

Fu and Kanc^"' considered a family of ID bulk-gapped 
Hamiltonians H{x) parametrized by a cyclic parameter 
t (i.e., H[t -|- r] = H[t]) subject to the constraint 

H[-t] = eH[t]0-\ (1) 

where 9 is the time-reversal operator. This can be un- 
derstood as an adiabatic pumping cycle, with t playing 
the role of time or pumping parameter. The constraint 
of Eq. (1) guarantees that the Hamiltonian H{x) is T- 
invariant at the points t = and t = T/2, while the 
T symmetry is broken at intermediate parameter values. 
If we also limit ourselves to Hamiltonians having unit 
period, so that H is invariant under x — > x + I, then 
the eigenstates may be represented by the periodic parts 
\unk) = e~^'^^\ipnk) of the Bloch states iV'nfc)- At t = 
and t = T/2 the Hamiltonian is time-reversal invariant 
and the eigenstates come in Kramers pairs, being degen- 
erate at fc = and k = n. 

Since the system is periodic in both k and t, the \unk) 
functions arc defined on a torus. Moreover, the system 
must also be physically invariant under a gauge transfor- 
mation of the form 

lUnk) = ^ U„in\Umk) (2) 
m 

where U{k,t) expresses the gauge freedom to 

choose Af representatives of the occupied space at each 



{k,t). We adopt a gauge that is continuous on the half- 
torus t e [0, T/2] and that respects T symmetry at t = 
and r/2 in the sense of Fu and Kane,^''' i.e., 

K-k) = -e'^'-'-Oluii}, 

\<'-k) = e'''-''0\ui,). (3) 

Here the occupied states n = 1, ...,7V have been relabeled 
in terms of pairs a = 1, A/'/2 and elements / and // 
within each pair. Note that Eq. (3) is a property which 
is not preserved by an arbitrary IA{N) transformation. It 
allows the Berry connection 

A{k) ^ i^{Unk\dk\Unk) (4) 
n 

to be decomposed as 

A{k)=A'{k) + A''{k) (5) 

where 

A''{k)=iY,{ui^\dk\ui^) (6) 

a 

and S = I, II. Having chosen a gauge that obeys these 
conventions at t = and r/2 and evolves smoothly for 
intermediate t,^^ the "partial polarizations"^'' 

P^ = ^fdkA'{k) (7) 

can be defined such that their sum is the total charge 
polarization^^ 

Pp^^jdkA{k) = pi + p;/. (8) 

Note that the total polarization is defined only modulo 
an integer (the quantum of polarization) under a general 
U{J\f) gauge transformation, while the "partial polariza- 
tion" is not gauge invariant at all. A quantity that is 
gauge-invariant is the change in total polarization during 
the cyclic adiabatic evolution of the Hamiltonian, and 
using Eq. (1) it follows that 

Pp(T)-P,(0) = C (9) 

where C is the first Chern number, an integer topolog- 
ical invariant corresponding to the number of electrons 
pumped through the system in one cycle of the pump- 
ing process. For a T-invariant pump that satisfies the 
conditions of Eq. (1), C must be zero. 

In order to describe the Z2 invariant of a T-symmctric 
system in a similar fashion, the "time reversal polariza- 
tion" was introduced as^^ 

Pe = Pp - P"- (10) 

Then the integer Z2 invariant can be written as 

A = Pe(r/2) - Pe(0) mod 2. (11) 
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To summarize, the Z2 invariant is well defined via 
Eq. (11) when the gauge respects T-symmetry at i = 
and T/2 and is continuous on the torus between these 
two parameter values. Note, however, that while such 
a gauge choice is possible on the half-torus even for the 
Z2-odd case (A=l), it can only be extended to cover the 
full torus continuously in the Z2-even case (A=0).^^'^*'^® 

B. Formulation in terms of Wannier charge centers 

Let us now rewrite Eq. (11) in terms of the Wannier 
charge centers (WCCs). By definition the Wannier func- 
tions (WFs) belonging to unit cell R are 

\Rn) = ^J dfce-*'=(«--)|u„fc). (12) 

The WCC Xn is defined as the expectation value x^ — 
{On\X\On) of the position operator x in the state |0n) 
corresponding to one of the WFs in the home unit cell 
R = 0. Equivalently,26>30 

dk(unk\dk\u„k)- (13) 

Except in the single-band case, the individual Xji are 
not independent of a general gauge transformation as in 
Eq. (2). However, the sum over all WCCs in the unit cell 
is a gauge-independet quantity (modulo a lattice vector, 
i.e., mod 1 in our notation). For the present purposes 
we adopt the gauge of Eq. (3) and construct WFs \Ra, S) 
by inserting litfj.) into the definition of Eq. (12). In this 
gauge 

xi = xi^ mod 1, (14) 

as follows from Eqs. (3) and (13) and use of the continuity 
condition Xa,-TT = Xa,7r + ^Trm, where m is an integer. 
Since we have also insisted on the gauge being continuous 
for t G [0,T/2], it is possible to follow the evolution of 
each WCC during the half-cycle. Taking into account 
that -£0.^ = (1/2^) /bz-^^ for S = I, II, Eq. (11) 
yields 

a a 

(15) 

Since the gauge is assumed to be smooth, the evolution 
of the charge centers must also be smooth. Being defined 
in this way, A is clearly a mod-2 quantity, and as shown 
in Ref. 24 it represents the desired Z2 invariant. 

However, if the gauge breaks T symmetry or it is not 
continuous in the half-cycle, Eq. (15) no longer defines 
a topological invariant. A discontinuity in the gauge in 
the process of the half cycle can change A by 1, so the 
mod-2 property is lost. Breaking T in the gauge choice 
means that the corresponding centers are not necessarily 
degenerate at t = and t — T/2. In fact, A can even 
take non-integer values in this case.^^ 



The above argument implies that in order to compute 
the Z2 invariant via Eq. (15), one needs a gauge that 
satisfies both T-invariance and continuity on the half- 
torus. We now argue that the gauge that corresponds to 
ID maximally localized WFs at each t has the desired 
properties, as long as these WFs are chosen to evolve 
smoothly as a function of t. The criterion introduced in 
Ref. 31 for constructing the maximally localized WFs was 
that the gauge choice should provide the minimum possi- 
ble quadratic spread = ^„[(0n|r^|0n) — (On|r|On)^]. In 
ID, the maximally localized WFs constructed according 
to this criterion are eigenstates of the position operator 
X in the band subspace.'^^''^^ Since this operator com- 
mutes with 9, its eigenvalues will be doubly degenerate 
and its eigenstates will come in Kramers pairs at i = 
and T/2. 

To prove continuity of this gauge in fc, let us briefly 
discuss how to enforce it on a fc-mesh fcj+i = kj + Ak by 
carrying out a multi-band parallel-transport construction 
along the Brillouin zone.'^^ At a given value of t, starting 

from fc = one constructs overlap matrices Mmn'^'^^^ = 
(wrnfcj Iwnfcj+i) m such a way that they are Hermitian. 
This can be done in a unique way by means of the singular 
value decomposition M = l/Eiy'', where E is positive 
real diagonal and V and W are unitary matrices. With 
this decomposition a unitary rotation of the states at 
fcj+i by WV'^ leaves M^'^^^'^J+i^ Hermitian. Repeating 
this procedure, one finds that states \ipnk) at fc = 27r are 
related to those at fc = by a unitary rotation A, whose 
eigenvalues An = e~"" give the ID maximally-localized 
WCCs Xn- The corresponding eigenvectors can be used 
to define a gauge that is continuous in fc for a given value 
of t. The continuity vs. t on the half-torus is achieved by 
tracing the evolution of the WCCs function of 

t, with the n'th state of the gauge constructed from the 
eigenvectors associated with the n'th smoothly evolving 

WCC Xn{t). 

Having established a particular gauge choice in which 
Eqs. (11) and (15) are valid, it is straightforward in prin- 
ciple to obtain the Z2 invariant. Indeed, Eq. (15) implies 
that the Z2 invariant can be determined simply by test- 
ing whether the WCCs change partners when tracked 
continuously from t=0 to T/2. This is the essence of our 
approach. We stress that no explicit construction of a 
smooth gauge on the half-torus is necessary; we simply 
track the evolution of the WCCs on the half-torus. 

In practice, when working on a discrete mesh of t values 
when many bands are present, it may not be entirely 
straightforward to enforce the continuity with respect to 
t. In the next section we present a simple and automatic 
numerical procedure that is robust in this respect, and 
use it to illustrate the calculation of the Z2 invariants for 
several materials of interest. 
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FIG. 1. Sketch of evolution of Wannier charge centers 
(WCCs) X vs. time t during an adiabatic pumping process. 
Regarding x £ [0, 1] as a unit circle and t £ [0,T/2] as a line 
segment, the cylindrical {x,t) manifold is represented via a 
sequence of circular cross sections at left, or as an unwrapped 
cylinder at right. Each red rhombus marks the middle of the 
largest gap between WCCs at given t. (a) Z2 insulator; WCCs 
wind around the cylinder, (b) Normal insulator; WCCs re- 
connect without wrapping the cylinder. 

III. NUMERICAL IMPLEMENTATION 

The method outlined above, in which the WCCs ob- 
tained with the ID maximally-localized gauge are used 
to compute the Z2 invariant via Eq. (15), can be im- 
plemented by plotting the WCCs at each point on the 
t mesh and then visually tracking the evolution of each 
WCC, as wc describe next in Sec. Ill A. However, we find 
that a more straightforward and more easily automated 
approach is to track the largest gap in the spectrum of 
WCCs instead. This gives rise to our proposed method, 
which is described in Sec. IIIB. 



A. Tracking WCC locations 

Let us first interpret Eq. (15) in terms of the winding 
of the WCCs around the BZ during the half-cycle t S 
[0,r/2]. Since the WCCs are defined modulo 1, one can 
imagine the x„ living on a circle of unit circumference, 
as illustrated in the left panels of Fig. 1. During the 
pumping process, the WCCs migrate along this circle. 
The system will be in the Z2-odd state ((5=1) if and only 
if the WCCs reconnect after the half cycle in such a way 
as to wrap the unit circle an odd number of times. 

Consider, for example, the case of only two occupied 
bands, as sketched in Fig. 1. The top panel shows the 
Z2-odd case; the blue and green arrows show the evo- 
lution of the first and second WCC from to (= 0) to ^4 
(= T/2), and they meet in such a way that the unit circle 



is wrapped exactly once. Correspondingly, as shown in 
the right-hand part of the figure, the WCCs "exchange 
partners" during the pumping process (i.e., two bands be- 
longing to the same Kramers pair at i = do not rejoin 
at t = For the Z2-even case shown in the bottom 

panel, by contrast, the unit circle is wrapped zero times, 
and no such exchange of partners occurs. 

If one has access to the continuous evolution of the 
WCCs vs. t, as shown by the solid blue and green curves 
in Fig. 1, this method works in principle for an arbitrary 
number of occupied bands (i.e., WFs per unit cell). An 
illustrative example with many bands appears in Fig. (1) 
of Ref. 33. Either the "bands" Xn exchange partners in 
going from t ~ to t = T/2 (0 = to = tt in their 
notation), or they do not, implying Z2 odd or even re- 
spectively. Equivalently, one can draw an arbitrary con- 
tinuous curve starting within a gap at t = and ending 
within a gap a.t t = T/2; the system is Z2-odd if this 
curve crosses the WCC bands an odd number of times, 
or Z2-even otherwise. 

In practice, however, one will typically have the WCC 
values only on a discrete mesh of t points, in which 
case the connectivity can be far from obvious. Certainly 
one cannot simply make the arbitrary branch cut choice 
Xn G [0, 1], sort the Xn in increasing order, and use the 
resulting indices to define the paths of the WCCs. This 
would, for example, give an incorrect evolution from ti to 
t2 in Fig. 1(b), since one WCC passes through the branch 
cut in this interval, apparently jumping discontinuously 
from the "top" to the "bottom" of the unwrapped cylin- 
der at right. (A similar jump happens again near ^3.) 

One possible approach is that of Ref. 33 mentioned 
above, i.e., to increase the t mesh density until, by vi- 
sual inspection, the connectivity becomes obvious. How- 
ever, this becomes prohibitively expensive in the first- 
principles context, since a calculation of many (typically 
10-30) bands would have to be done on an extremely fine 
mesh of t points. It is typical for some of the WCCs 
to cluster rather closely together during part of the evo- 
lution in t; if this clustering happens near the artificial 
branch cut, it can become very difficult to determine the 
connectivity from one t to the next, even if a rather dense 
mesh of t values is used. Moreover, an algorithm of this 
kind is difficult to automate. For these reasons, we find 
that the direct approach of plotting the evolution of the 
WCCs is not a very satisfactory algorithm for obtain- 
ing the topological indices, at least in the case of a large 
number of occupied bands. 



B. Tracking gaps in the WCC spectrum 

Here we propose a simple procedure that overcomes 
the above obstacles, allowing the Z2 invariant to be com- 
puted in a straightforward fashion. The main idea is to 
concentrate on the largest gap between WCCs, instead of 
on the individual WCCs themselves. As explained above 
and illustrated by the red dashed curve in Fig. 1, the 



path following the largest gap in a;„ values (with vertical 
excursions at critical values of t) crosses the a;„ bands a 
number of times that is equal, mod 2, to the Z2 invariant. 
Our approach, in which we choose this path as an espe- 
cially suitable one for discretizing, can be implemented 
without reference to any branch cut in the determination 
of the Xn, allowing the Z2 invariant to be determined 
from the flow of WCCs on the cylindrical (x, t) manifold 
directly. 

As in Fig. 1, we again consider a set of M circular sec- 
tions of the cylinder that correspond to the pumping pa- 
rameter values = T{m - 1)/2M, where m e [0, M]. 
At each tm we define z*^™^ to be the center of the largest 
gap between two adjacent WCCs on the circle. (If two 
gaps are of equal size, either can be chosen arbitrarily.) 
For dcfinitcncss we choose z(™) e [0, 1), but as we shall 
sec shortly, the branch choice is immaterial. In the con- 
tinuous limit 7\/ — > 00, z(f) takes the form of a series of 
path segments on the surface of the cylinder, with dis- 
continuous jumps in the x direction at certain critical pa- 
rameter values tj . Our algorithm consists in counting the 
number of WCCs jumped over at each tj, and summing 
them all mod 2. As becomes clear from an inspection 
of Fig. 1 and similar examples of increasing complexity, 
the WCCs exchange partners during the evolution from 
t=0 to T/2 only if this sum is odd, so that this sum 
determines the Z2 invariant of the system. 

The approach generalizes easily to the case of discrete 



z^"'\ Let A„, be the number of WCCs x^ 



(m+l) 



that ap- 



pear between gap centers z'™' and 2;('"+i), mod 2. As we 
shall see below, this can be computed in a manner that is 
independent of the branch cut choices used to determine 
the x™ and z^™). Then the overall Z2 invariant is just 



A 



M 

E 

m=0 



A,„ mod 2. 



(16) 



This argument is illustrated in the right-hand panels 
of Fig. 1 for the two band-case and M = 4. The rect- 
angles represent the surface of the cylinder in the pa- 
rameter space, and should be regarded as glued along 
the longer sides. The circles correspond to i^™-* values, 
while each red rhombus represents the center z'™' of the 
largest gap between il™'' values. In Fig. 1(a) there is 
one jump that occurs between m=2 and m=3, in which 
one WCC is jumped over; thus, A,„ = except for 
A2 = 1, giving A=:l. In Fig. 1(b), on the other hand, 
there are two jumps, once between m—1 and m=2 and 
again between m—2 and m=3, so that Ai = A2 = 1 and 
A = (mod 2). 

We now show how the A™ can be computed straight- 
forwardly in a manner that is insensitive to the branch- 
cut choices made in determining the and z^'^\ We 
use the fact that the directed area of a triangle defined 
by angles 0i, (1)2, and 03 on the unit circle is"^^ 

02,03) sin(02 -0i) + sin(03-02) + sin(0i -03). 

(17) 




FIG. 2. Sketch illustrating the method used to determine 
whether ii""^^' lies between z'™' and 2'™+^' in the counter- 
clockwise sense when mapped onto the complex unit circle, 
(a) Yes, since the directed area of the triangle is positive, (b) 
No, since it is negative. 



Therefore the sign of .g(0i, 02, 03) tells us whether or not 
03 lies "between" 0i and 02 m the sense of counterclock- 
wise rotation. Identifying 0i = 2'kz^"^\ 02 = 27rz(™+-'^^ 
and 03 = 2iTx^^^\ as in Fig. 2, it follows that 



(-1)^ 



n ^g^^ 



n=l 



g(27^z('"^27rz^™+^^27r4' 



(m+l) 



(18) 

where sgn(a;) is the sign function. The A„i defined in 
this way is precisely the needed count of WCCs jumped 
over, mod 2, in evolving from m to to -|- 1. 

As a last detail, we discuss the case of possible de- 
generacies between the three arguments of .g(0i, 02, 03). 
First, note that 2:(™+i) = xll^'^^'' is impossible, since 
^{m+i) jg (definition in a gap between Xn^'^'^^ values. If 
the mesh spacing in t is fine enough, then by continuity 
we expect that z^™-' = xl'"^^'' will also be unlikely. It is 
recommended to test whether these values ever approach 
within a threshold distance, and restart the algorithm 
with a finer t mesh if such a case is encountered; two 
cases of this kind are discussed later in Sec. IV. Finally, 
it can happen that zf^™^ = z^'^'^^\ In this case, the 
signum function (which technically assigns value to ar- 
gument 0) should be replaced in Eq. (18) by a function 
that returns s whenever z^™^ = 2:(™+^\ where s is chosen 
once and for all to be either -1-1 or —1. Since the same 
degeneracy appears in every term of the product over TV 
factors in Eq. (18), where Af is even, the choice of s is 
arbitrary as long as it is applied consistently. 

The above-described algorithm, based on Eqs. (16-18), 
constitutes one of the principal results of the present 
work. The implementation of this algorithm is straight- 
forward, and allows for an efficient and robust determi- 
nation of the Z2 invariant even when many bands are 
present, and even for only moderately fine mesh spac- 
ings. In Sec. IV, we will demonstrate the successful ap- 
plication of this approach to the calculation of the strong 
and weak topological indices of some real materials. 
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C. Application to 2D and 3D T-invariant insulators 

As pointed out in Rcf. (24). the pumping process dis- 
cussed above for a ID system is the direct analogue of 
a 2D T-invariant insulator, i.e., one whose Hamiltonian 
is subject to the condition k) = 9^^H{k)9. To see 
this, let k = ^jfcibi/27r, where bi and b2 have been 
chosen as primitive reciprocal lattice vectors. Then we 
can let ki and fc2 play the roles of k and t respectively. 
Just as H{k, t) displays T symmetry of H{x) at t=0 and 
T/2, so H{ki,k2), regarded as the Hamiltonian H{xi) 
of a fictitious ID system for given A;2, is 7~- invariant at 
k2 = and tt. The Wannier functions of the effective 
ID system can be understood as "hybrid Wannier func- 
tions" that have been Fourier transformed from k space 
to r space only in direction 1, while remaining extended 
in direction 2. The topological Z2 invariant of the 2D 
system can therefore be determined straightforwardly by 
applying the approach outlined above. 

A topological phase of a 3D T-symmetric insulator is 
described by one strong topological index and three 
weak indices vi, 1/2, and ly^.^-'^-^^ These indices may be un- 
derstood as follows. Again letting k = fcib,;/27r, there 
are eight 7~-invariant points ^(ni,7i2,n3)j where n,; = or 1 
denotes ki = or tt respectively. These eight points may 
be thought of as the vertices of a parallelepiped in recip- 
rocal space whose six faces are labeled by ni=0, n2=0, 
713=0, ni=l, 712=1, and 713=1. On any one of these six 
faces, the Hamiltonian -ff(k), regarded as a function of 
two k variables, can be thought of as the Hamiltonian of 
a fictitious 2D T-symmetric system, and the argument 
of the previous paragraph can thus be applied to each of 
these six faces separately. The three weak indices Vi=i.2,z 
are defined to be the Z2 invariants associated with the 
three surfaces 77i=l, 772=1, and 773=1."^ These weak in- 
dices obviously depend on the choice of reciprocal lattice 
vectors. The strong index vq is the sum (mod 2) of the 
I12 invariants of the 77^=0 and nj=l faces for any one of 
the j (implying some redundancy among the six indices) ; 
it is also a Z2 quantity, but is independent of the choice 
of reciprocal lattice vectors. ^'^ 

Thus, a complete topological classification in 3D, given 
by the index lyo', {i'ih'21^3) , can be obtained by applying 
our analysis to each of these six faces in the 3D Brillouin 
zone. Note that in general, this determines the strong 
index with some redundancy, providing a check on the 
internal consistency of the method. However, symmetry 
considerations often play a role. For systems having a 
3-fold symmetry axis, for example, one typically needs 
to compute the Z2 index on only two faces, as we shall 
see below. 



IV. APPLICATION TO REAL MATERIALS 

In this section we discuss the application of the above- 
described method to real materials. First, we illustrate 
the validity of the approach for centrosymmetric Bi and 



212803, where weak and strong indices may alternatively 
be computed directly from the parities of the occupied 
Kramers pairs at the eight T-invariant momenta.^ We 
then apply the method to noncentrosymmetric crystals 
of GeTe and strained HgTe, showing that the first is a 
trivial insulator, while the latter is a strong topological 
insulator under both positive and negative strains along 
[001] and under positive strain along [111]. 

The calculations were carried out in the frame- 
work of density-functional theory"^^ using the local- 
density approximation with the exchange and corre- 
lation parametrized as in Ref. 37. We used HGH 
pscudopotentials'^* with semicore Sd-states included for 
Hg, while for all other elements only the s and p va- 
lence electrons were explicitly included. The calculations 
were carried out using the ABINIT code package"^^'^^ 
with a 10 X 10 X 10 k-mesh for the self-consistent 
field calculations and a 140 Ry planewave cutoff. The 
spin-orbit interaction was included in the calculation 
via the HGH pseudopotentials. Note that the over- 
lap matrices Mmn^'*^^ defined in Sec. II B, are the 
same as those needed for the calculation of the elec- 
tric polarization^^ or the construction of maximally- 
localized Wannier functions, and arc thus readily avail- 
able in many standard ah initio code packages including 
ABINIT. 



A. Centrosymmetric materials 

We start by illustrating the method with the examples 
of Bi and Bi2Se3. Although Bi is a semimetal, its ten 
lowest-lying valence bands are separated from higher ones 
by an energy gap everywhere in the BZ, so in this case the 
topological indices describe the topological character of a 
particular group of bands. Since this is not the occupied 
subspace of an insulator, these topological indices are not 
"physical," but it is still of interest to compute them and 
compare with methods based on the parity eigenvalues.® 
According to the latter approach, the group of ten lowest- 
lying bands of Bi was shown to be topologically trivial.® 
Bi2Se3, on the other hand, is a true insulator, and the 
parity approach demonstrated that it is a strong topo- 
logical insulator.'' 

Bi and Bi2Sc3 both belong to the rhombohedral space 
group Rim (#166), which has a 3-fold rotational axis. 
Thus, it is enough to compute only one weak Z2 index, 
say for 711 = 1, since all three of them are equal by sym- 
metry. To get the strong index, one just needs to compute 
just one more of the Z2 invariants, say for 771 =0. 

Our results for Bi, obtained with the lattice parameters 
used in previous studies,"*^ arc presented in Fig. 3. Panels 
(a) and (b) show the determination of the Z2 invariant 
at 77i=0 and rii = l respectively, with ^2 treated as the 
pumping parameter (like t) for an effective ID system 
with wavevector k-^. The k2 axis was initially discretized 
into ten equal intervals (777 = 1, 10) running from to 
TT. but for reasons discussed below an extra point (number 



7 




0.5 

(b) 



0.0 



-0.5 







o 


o 


o 


u 


u 


o 


o 




▼ 










o 





o 


o 


o 


o 


o 


o 

§ 


8 




o 







o 





o 


o 


o 


o 







o 


♦ 






o 


o 


o 


o 


o 


o 


o 





o 


o 


o 


o 






o 


o 


o 


o 


o 


o 


o 


o 


o 


o 






♦ 


♦ 


♦ 


♦ 


♦ 


♦ 


♦ 


♦ 












o 


o 
o 






o 
o 


o 
o 


o 
o 


oo 


oo 


o 
o 


o 
o 


o 
o 


o 





















♦ 






♦ 




o 


o 


o 


o 


o 


o 







o 










o 





o 


o 


o 


o 





§ 

o 


8 




o 






o 











o 













1 




1 

1 


1 

2 


3 


4 


5 


1 

6 


1 

7 


1 

8 


9 


1 

10 


T 

11 




8 9 10 



FIG. 3. Evolution of Bi WCCs x„ (circles) in the rs direction 
vs. k2 at (a) fci=0; (b) ki—ir. Red rhombus marks midpoint 
of largest gap. k2 is sampled in ten equal increments from 
to TT, except that an extra point is inserted midway in the last 
segment in panel (b) (see text). 

10 on the horizontal axis of the plot) was inserted midway 
in the last segment to make a total of eleven m values 
in Panel (b). As noted above, we are treating a group of 
ten valence bands labeled by n, so we have an array of 
WCC values il™'' whose values are indicated by the black 
circles in the plot. These form Kramers pairs at k2—0 
and TT, but not elsewhere. Each red rhombus indicates 
the center z*^™) of the largest gap between adjacent xi™'' 
values, as discussed in Sec. III. 

Looking first at Fig. 3(a), we see that the gap center 
jumps over one WCC at m=l, and then over three WCCs 
at m=7, for a total of four, which is even. In Fig. 3(b) 
we get a total of 2 + 7 + 3 + 4= 16 jumps, which is again 
even. The visual determinations of the number of jumped 
bands is confirmed by the application of the automated 
procedure of Eqs. (16-18). Thus, both Z2 indices are 0, 
and the 3D index is 0; (000), indicating a normal band 
topology as anticipated. 

We now discuss the above-mentioned insertion of one 
extra k2 point in Fig. 3(b). This was necessary because 
the gap center z^^^ at k2 = O.Qtt had almost the same 
value as one of the WCC values a.t k2 — tt (now labeled as 
'11' on the horizontal axis), making it ambiguous whether 
or not that a;„ value should be counted as one of the ones 
that has been jumped over. To resolve this difficulty, we 
included an extra step at fc2 — 0.957r (now labeled as '10' 
on the horizontal axis). The reason for the fast motion 
of the WCC in this case is that the minimum gap to the 
next higher (eleventh) band becomes rather small near 

k2 = TT. 

Note that the detection of this kind of problem does 



FIG. 4. Evolution of Bi2Se3 WCCs x„ (circles) in the ra 
direction vs. ^2 at (a) fci=0; (b) ki—ir. Red rhombus marks 
midpoint of largest gap. ^2 is sampled in ten equal increments 
from to TT. 

not have to be done by visual inspection, but can be au- 
tomated in the context of Eqs. (16-18). As already men- 
tioned in Sec. IIIB. we simply test whether any ii™^^' 
approaches within a certain threshold of z*^™^ (mod 1); 
if so, we flag the interval in question for replacement by 
a finer mesh. Still, it is recommended to choose a mesh 
that is fine enough so that this threshold is rarely encoun- 
tered, with a finer mesh recommended in cases where the 
minimum band gap is small. 

The analysis of the same ni=0 and ni=l faces for the 
28 WCCs of BiaSes is illustrated in Fig. 4. The experi- 
mental lattice parameters''^ were used. Here there are no 
jumps over WCCs except for a single one in the very first 
step in the top panel (ni=0). It follows that the topo- 
logical index is 1; (000), in accord with previous studies.'' 

B. Noncentrosymmetric materials 

We now proceed to systems without inversion symme- 
try, which are the principal targets of our method since 
an analysis based on parity eigenvalues is not possible. 

GeTe belongs to the rhombohedral R3m space group 
(#160) and has no inversion symmetry, although like 
Bi and Bi2Se3 it has a 3-fold rotational symmetry, so 
that only two reciprocal-space faces have to be studied. 
The experimental lattice parameters''* were used, and the 
evolution of the 10 WCCs is presented in Fig. 5 follow- 
ing similar conventions as for Bi and Bi2Se3. For both 
faces Eq. (18) gives a trivial Z2 index, with the center 
of the largest gap making no jumps, so that GeTe is in 
the topologically trivial state 0; (000). This result could 
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FIG. 5. Evolution of GeTe WCCs S„ (circles) in the ra di- 
rection vs. k2 at (a) fci=0; (b) ki—n. Red rhombus marks 
midpoint of largest gap. k2 is sampled in ten equal incre- 
ments from to TT. 



have been anticipated from the fact that the spin-orbit 
interaction in GeTe is weak, as reflected in the approxi- 
mate pairwise degeneracy of the WCCs throughout the 
evolution. 

Finally, let us consider the more interesting case of 
strained HgTc. In the absence of strain this is a zero- 
band-gap material. Any anisotropic strain breaks the 
four- fold symmetry at F, making it possible that the gap 
might open. Based on an adiabatic continuity argument, 
HgTe was predicted to be a strong topological insula- 
tor under compressive strain in the [001] direction.^ This 
was later verified with tight-binding calculations.^'*^ Ap- 
plication of our approach to HgTe under uniaxial strain 
also confirms that HgTc is a strong topological insulator, 
with index 1; (000), under both positive and negative^ 
2% strains along the [001] direction (not shown). This 
means that although the positive-strain and negative- 
strain states are separated by a gap closure at zero strain, 
there is no topological phase transition associated with 
this gap closure. 

We also studied strains in the [111] direction. Un- 
der compressive strains of —2% and —5% the system 
becomes metallic and the direct band gap vanishes, so 
that no topological index can be associated with the oc- 
cupied space. Under tensile strain of 4-2% we find that 
HgTc becomes a narrow-gap semiconductor with an in- 
direct energy gap of Eg = 0.054 eV, while for -1-5% strain 
it becomes metallic. Even at -1-5% strain, however, the 
lowest 18 bands remain separated from higher ones by an 
energy gap at all k, so that, as for Bi, one can still assign 
a topological index to this isolated group of bands. The 
computed band structures for both cases are illustrated 
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FIG. 6. Band structure along high-symmetry lines of the 
undistorted FCC structure for HgTe under tensile strain in 
the [111] direction, (a) +2% strain, (b) -1-5% strain. 
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FIG. 7. Evolution of WCCs for HgTe under +2% strain in 
the [111] direction. WCCs x„ (circles) in the rs direction are 
plotted vs. k2 at (a) fci=0; (b) ki—n. Red rhombus marks 
midpoint of largest gap. k2 is sampled in ten equal increments 
from to TT, except that an extra point is inserted midway in 
the first segment in panel (a) (see text). 



in Fig. 6 along lines connecting the high symmetry points 
of the undistorted FCC structure. 

The space group of [lll]-strained HgTe is rhombohe- 
dral R3m (#166), the same as for GeTe, so that again 
only two Z2 indices need to be calculated. The results of 
our WCC analysis for the case of -1-2% strain arc shown 
in Fig. 7. We find Z2=l and Z2=0 for ni=0 and 722=1 
respectively, so that the topological class is 1; (000). The 
behavior in Panel (b) is rather uninteresting, since the 
gap is large everywhere on the n2=l face. However, in 
Panel (a) we again find an example of a rapid change 
of WCCs with fc2, which was repaired by inserting an 
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extra point (the one now labeled '1' on the horizontal 
axis) at k2 = O.OStt. Actually, we anticipated the need 
for this denser sampling for small ^2 from the fact that 
the zero-strain gap closure occurs at F, so that a delicate 
dependence on k near the BZ center was expected. 

V. SUMMARY AND CONCLUSIONS 

We have proposed a new approach for calculating topo- 
logical invariants in T-invariant systems. The method 
is based on following the evolution of hybrid Wannier 
charge centers, and is very general, being easily applica- 
ble in both tight-binding and DFT contexts. The needed 
ingredients are the same as those needed for the cal- 
culation of the electric polarization or the construction 
of maximally-localized Wannier functions, and are thus 
readily available in standard code packages. The present 
algorithm is relatively inexpensive, however, because the 



analysis is confined to a small number of 2D slices of 
the 3D Brillouin zone. The method is easily automated 
and remains robust even when many bands are present. 
We hope that our method can help to make the search 
for topological phases in noncentrosymmetric materials a 
routine task, and that it will lead to further progress in 
this rapidly developing field. 

Note: In the final stages of preparing this manuscript, 
we became aware of independent work by Yu et al.^^ that 
is closely related. These authors carry out a similar anal- 
ysis based on WCCs, but without the automated analysis 
described in our Sec. III. 
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